Cluster Abundance in f(R) Gravity Models 
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As one of the most powerful probes of cosmological structure formation, the abundance of massive 
galaxy clusters is a sensitive probe of modifications to gravity on cosmological scales. In this paper, 
we present results from ./V-body simulations of a general class of f(R) models, which self-consistently 
solve the non-linear field equation for the enhanced forces. Within this class we vary the amplitude 
of the field, which controls the range of the enhanced gravitational forces, both at the present epoch 
and as a function of redshift. Most models in the literature can be mapped onto the parameter 
space of this class. Focusing on the abundance of massive dark matter halos, we compare the 
simulation results to a simple spherical collapse model. Current constraints lie in the large-field 
regime, where the chameleon mechanism is not important. In this regime, the spherical collapse 
model works equally well for a wide range of models and can serve as a model-independent tool 
for placing constraints on f(R) gravity from cluster abundance. Using these results, we show how 
constraints from the observed local abundance of X-ray clusters on a specific f(R) model can be 
mapped onto other members of this general class of models. 
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I. INTRODUCTION 

The abundance of massive galaxy clusters provides a 
unique test of gravity on cosmological scales [l|-|3[ ■ Once 
constrained to expansion history data, modified gravity 
explanations of the cosmic acceleration generically pre- 
dict very different effects on the growth of cosmological 
structure than spatially smooth dark energy like the cos- 
mological constant. Moreover as highly non-linear ob- 
jects, clusters provide a testing ground for the non-linear 
interactions of viable theories where gravity becomes in- 
distinguishable from General Relativity locally. 

In the so-called f(R) class of models (see [U, Q and 
references therein) the modification to gravity arises from 
replacing the Einstcin-Hilbert action by a function of the 
Ricci or curvature scalar R 6 8]. These models possess 
an extra scalar degree of freedom fn = df/dR which 
mediates a 4/3 enhancement of gravitational forces on 
scales below the Compton wavelength or range associated 
with its mass. 

This enhancement changes the abundance of rare dark 
matter halos associated with clusters of galaxies. Mea- 
surements of the cluster abundance provide the current 
best cosmological constraints on f(R) models [HQ. On 
the other hand, in order to hide these enhancements from 
local tests of gravity, viable f(R) models employ the 
chameleon mechanism which allows the Compton wave- 
length to shrink in regions with deep gravitational po- 
tential wells 0, [l(j ■ Cosmological simulations including 
the chameleon effect are required to explore the impact 
of these modified forces on the cluster abundance. These 
have so far been performed for only a specific form of 

f(R) [nun. 

In fact, the relationship between the Compton wave- 
length, chameleon threshold and their respective evolu- 
tion with redshift depends on the functional form of f(R). 
In this paper, we explore the dependence of the cluster 



abundance on the functional form of f(R) in order to 
place more robust constraints on the whole class of mod- 
els. 

In ijnl we review the phenomenology of f(R) models, 
simulation technique and spherical collapse modeling as 
well as show that a general class of broken power law 
models introduced in Ref. [3] covers most cases of cos- 
mological interest. In §1111 we study the enhancement of 
the cluster abundance in these models and obtain con- 
straints from the local X-ray sample. We discuss these 
results in EIIVI 



II. METHODOLOGY 



We begin in i jll Al with a review of f(R) models. In 
Bl we discuss the numerical TV-body simulations from 
which wc extract the cluster abundance enhancements. 
In mi Bl wc discuss the semi-analytic modeling of these 
results with spherical collapse collapse calculations. 



A. Models 

In the f(R) model, the Einstein-Hilbert action is aug- 
mented with a general function of the scalar curvature 
R, 
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Here and throughout c = h = 1. Gravitational force 
enhancements are associated with an additional scalar 
degree of freedom, the chameleon field = df/dR, 
and have a range given by the comoving Compton wave- 
length Ac = a~ x {idf r/ dR) 1 / 2 . This additional attrac- 
tive force leads to the enhancement in the abundance of 
rare massive dark matter halos described below. The 
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FIG. 1: Rcdshift evolution of the chameleon field fn (top panel) 
and Compton wavelength (bottom panel) in the background 
for the broken power law class of models. As the scaling index 
n increases the field amplitude becomes increasingly suppressed 
leading to stronger chameleon effects for the same gravitational 
potentials of clusters. The Compton wavelength for a fixed field 
amplitude today remains relatively constant at z < 1 and then also 
becomes increasingly suppressed with n. An alternative class of 
models specified by the expansion history and Compton wavelength 
Bq parameter is also shown for comparison (dashed lines). 



second important property of such models is the non- 
linear chameleon effect which shuts down the enhanced 
forces in regions with deep gravitational potential wells 
compared with the field at the background curvature R, 

1*1 > \fn{R)\- 

Given that different models for f(R) produce differ- 
ent scalings of the Compton wavelength and chameleon 
threshold with curvature and hence implicitly with red- 
shift and the degree of non-linearity, we wish to explore 
the dependence of the halo abundance with variations in 
the form of f(R). 

We therefore choose a class of models where the sca ling 
index with curvature can vary as a broken power law [14| 
such that 
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with two free parameters, A, /x 2 for each value of the 
scaling index n. Note that as R — > 0, f(R) —> 0, and 
hence these models do not contain a cosmological con- 
stant. Nonetheless as R ^> fi 2 , the function f(R) can be 
approximated as 



f(R) 
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with /ro = — 2nA/x 2 "/_R l+1 replacing /x as the second pa- 
rameter of the model. Here we define -Ro = R(z = 0), so 
that fuo = ffi(Ro), where overbars denote the quantities 
of the background spacetime. Note that if \/ro\ *C 1 the 
curvature scales set by A = O(i?o) and differ widely 
and hence the R 3> /x 2 approximation is valid today and 
for all times in the past. 

The background expansion history mimics ACDM with 
A as a true cosmological constant to order f R0 . Therefore 
in the limit \ f R0 \ < 10" 2 , the f(R) model and ACDM are 
essentially indistinguishable with geometric tests. On the 
other hand, the field amplitude parameter (fm) controls 
the range of the force modification and the chameleon 
mechanism. With the functional form of Eq. ([3]), the 
comoving Compton wavelength becomes 
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with a value at the background curvature today Rq = 
3Hg(4-3£l m ) of 



A 



CO 



16.6a 



'I/. 



7?0 



10" 



3a 



-/i _1 Mpc, 



(5) 



assuming a flat universe. As the scaling index n increases, 
the Compton wavelength today increases given the same 
background field amplitude today ffj . Conversely as n 
increases, force modifications at high rcdshift versus to- 
day decrease and the chameleon mechanism extends to 
shallower potential wells. Thus the net effect is a fairly 
weak dependence of Ac on n at z < 1. In Fig. [TJ we 
show the evolution of the background field and Comp- 
ton wavelength for the ft m = 1 — Q,\ = 0.24, h = 0.73 
cosmology that we simulate below. 

This set of broken power law models covers the cos- 
mological phenomenology of most viable f(R) models. 
For example the models of Ref. [l5| compose a subset of 
this class. It also has sufficient flexibility to bracket the 
behavior of models where the combination of a specific 
expansion history [T|| [TtJ and the Compton wavelength 
today fixes the form of f(R) [3- For the ACDM expan- 
sion history and a dimensionless Compton wavelength 
parameter 



_ df R /dR H 



2.1ft 



-0.76 



/ho | 



(6) 



where ' = d/dlna, the redshift evolution goes from n ~ 4 
at low curvature and redshift to n ~ 0.13 at high cur- 
vature and redshift (see Fig. [1]). For a fixed \fno\ the 
amount of linear growth at z = in the Bq model is 
smaller than in the n — 1 model and this must be borne 
in mind when c omp aring constraints between the two 
models (cf. [3, H, [T^). 

Likewise these models have stronger chameleon effects 
at z < 1 than the n = 1 broken power law model. A 
similar caveat applies to models with exponential rather 
than power law suppression of the field with curvature 
(e.g. E3). 
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B. Simulations 



TABLE I: Summary of simulations used for this work 



We conduct TV-body simulations of these broken power 
law models with a particle- mesh relaxation code [IJ, LL2| ■ 
Briefly, at each time step we first solve the non-linear 
field equation for the field fluctuation, 



V 2 % = - [SR(f R ) - ^GS Pra ] , 



(7) 



using a multigrid relaxation scheme. Here coordinates 
are comoving, Sf R = Jr{R) - fn(R), SR = R - R, 
Spm = Pm — Pm- The Sfn field fluctuation then acts 
as an additional source to the gravitational potential, 



V 2 * = 4nGa 2 Sp„ 



1. 



-y 2 sf R . 
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This linear equation for ^ is solved via a fast Fourier 
transform. Once ^ is known on the mesh, particles are 
moved in the usual way. 

Since the field equation implies that spatial variations 
in dfa will be of order the gravitational potential, there 
are two regimes of interest. In the large- field regime, 
the background value fn(R) is large compared with the 
gravitational potentials of structure, and the field equa- 
tion ([7]) can be linearized via 
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dR 
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where Ac (a) is evaluated at the background curvature 
R(a). In this case the joint solution of the Poisson and 
field equations in Fourier space is 



fc 2 * 
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Hence the background Compton wavelength sets the 
global range of the enhanced gravitational force. We call 
this the no chameleon case and for comparison conduct 
separate simulations employing Eq. (I10|) . 

In the small-field regime, /r(-R) is comparable to or 
smaller than typical gravitational potentials of structure, 
so that the curvature changes non-linearly with the field. 
In other words the Compton wavelength depends on the 
local curvature or field Ac = Ac(a,x). Field fluctuations 
saturate in deep gravitational potential wells (/r — > 0), 
leading to an equilibrium solution SR = 8irGSp m and a 
suppression of non-Newtonian forces. 

We use simulations of three different box sizes 
(400, 256, 128Mpc//i), and 6 simulations for each box size 
and model. The runs and models as well as mass resolu- 
tion for each box are summarized in Table HI To reduce 
the effect of sample variance, we compare each f(R) sim- 
ulation run to a ACDM simulation with the same initial 
conditions, i.e. the same initial density field drawn from 



an initial power spectrum with A s 
k = 0.05MPC 1 and n s = 0.958. 



(4.73 x 10- 5 ) 2 at 





|/ho| 


400 


256 


Mpc) 
128 


#full 


1(T 4 (n=l, 2) 
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runs 


3 • 1CT 6 (n=2) 


6 


6 


6 




1(T 6 (n=l) 
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#no 


10" 4 (n=l/2, 1, 2, 4, 8) 


6 


6 


6 


cham. runs 


3 ■ 1CT 6 (n=2) 


6 


6 


6 




1(T 6 (n=l) 


6 


6 


6 


ACDM 





6 


6 


6 


M h , n 


lin (10 12 /i- 1 M Q ) 


204 


53.7 


6.61 



We measure the mass function from the simulations 
using the methodology described in [l3| and refer the 
reader to details therein. Briefly, we identify halos using 
a spherical overdensity criterion of A = 200 with respect 
to the mean density and quantify the mass function en- 
hancements of the f(R) models over ACDM with the 
same initial conditions. To reduce the effect of shot noise 
we bin results into coarse mass intervals corresponding 
to approximately an e-fold (A In A/200 = 1-04). Further- 
more, due to resolution effects, we only utilize halos that 
contain at least 800 particles corresponding to the mini- 
mum mass given in Tab. [I] 

We estimate sampling errors via bootstrap resampling. 
Note that due to our limited number of realizations, these 
errors might be underestimated at high masses where 
halos are rare and fluctuations are significant. 



C. Spherical Collapse Predictions 

Since the large field regime is where the current local 
cluster abundance measurements constrain f(R) models 
[HHf , characterizing this regime in a way that does not re- 
quire simulations of each model is important. We briefly 
review a method utilizing spherical collapse introduced 
in Ref. [3 

The Sheth-Tormen description for the comoving num- 
ber density of halos per logarithmic interval in the virial 
mass M v is given by 



dn 
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EnLfty)-*!- 

M v V 'd\nM v 



11 



where the peak threshold v = S c /a(M v ) and 



vj{v) = A^au 2 [l + K 2 )" p ] exp[-a^ 2 /2] . (12) 

Here o~(M) is the variance of the linear density field con- 
volved with a top hat of radius r that encloses M = 
47rr 3 p m /3 at the background density 



a 2 (r) 



d 3 k 



\W(kr)\ 2 P h (k), 
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where -Pl(&) is the linear power spectrum and W is the 
Fourier transform of the top hat window. The normal- 
ization constant A is chosen such that / duf(v) = 1. The 
parameter values of p = 0.3, a = 0.75, and S c = 1.673 
for the spherical collapse threshold have previously been 
shown to match simulations of ACDM at the 10 — 20% 
level. The virial mass is defined as the mass enclosed 
at the virial radius r v , at which the average density is 
A v times the mean density. The virial mass can then be 
transformed to alternate overdensity criteria assuming a 
Navarro-Frenk- White density profile pl| . 

Spherical collapse can also provide a model for the 
mass function enhancement measured in the f(R) N- 
body simulations [T^j. The mass function calculation 
again uses the Sheth-Tormen form of Eq. pT|) but with 
the linear power spectrum for the f(R) model in Eq. (|13|). 
and two limiting cases for the spherical collapse param- 
eters. In one case, we simply assume that the spher- 
ical perturbation considered is always larger than the- 
Compton wavelength of the Jr field, so that gravity is 
GR throughout, and the spherical collapse parameters 
are unchanged. In the second case, we assume that the 
perturbation is always smaller than the local Compton 
wavelength in spite of the redshift evolution of the back- 
ground Compton wavelength and chameleon mechanism 
(see Fig.[T]). Hence forces arc simply universally enhanced 
by 4/3. In both cases, we use the modified linear force 
calculation for the linear power spectrum and o~(M) via 
Eq. (|11[) . Hence, unmodified spherical collapse param- 
eters does not equate to unmodified spherical collapse 
predictions. 

The values of the resulting linear collapse threshold S c 
and virial overdensity A v are summarized in Table [TTJ We 
use the GR values to calculate the mass function Eq. (fTTj) 
in terms of virial mass M v (A/ v = M A J for ACDM, 
and correspondingly for f(R) with either set of collapse 
parameters. We then rescale both mass functions to our 
adopted mass definition M200 and convolve them with 
the mass binning used in the simulations before taking 
the ratio. 



III. CLUSTER ABUNDANCE 

With the f(R) simulations described in Tab. |U we 
can now test the model dependence of the cluster abun- 
dance enhancement as well as the accuracy of the model- 
independent spherical collapse technique described in the 
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FIG. 2: Mass function enhancement at z = in large field 
I/-R0I = 10~ 4 models for scaling index of n = 1,2 relative to 
ACDM. Here and in the following figures, the no-chameleon re- 
sults have been displaced horizontally for clarity. Enhancement 
depends mainly on mass due to the increasing rarity of high mass 
halos. As n increases, the enhancement drops only moderately 
given the small change in the background Compton wavelength at 
z < 1, consistent with only a small contribution from the non- 
linear chameleon effect. The spherical collapse predictions (shaded 
range) capture these qualitative trends and provide conservative 
lower limits to the enhancement. 

previous section. In mil Al we discuss the large field 
regime relevant for current constraints from clusters. In 
Wl Bl wc evaluate the impact of the non-linear chameleon 
mechanism in the small field regime. Finally we show 
how constraints on one f(R) model can be transformed 
to another using simulation calibrated spherical collapse 
methods in ^HICl 



A. Large Field Regime 

In Fig. [SJ we show the mass function enhancements for 
a large field case \fno\ = 10 -4 for n = 1,2. Note that 
we plot the data points at the center of each mass bin, 
while the average mass of halos within the bin is generally 
smaller than that due to the steepness of the mass func- 
tion. The spherical collapse predictions are convolved 
with the mass bin and hence take into account this ef- 
fect. The uppermost mass bin extends to infinite mass 
so as to include all remaining halos but is still plotted at 
AlnM2oo = 1-04 above the previous bin. 

As the mass increases and halos become rarer in the 
ACDM simulations, the fractional impact of the force en- 
hancement on cluster abundance increases. Relative to 
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FIG. 3: Mass function enhancement at z = as a function of 
scaling index n for the mass bin centered at M200 = 5.3x 10 14 Mq//i 
for large field models |/ijol = 10~ 4 . Spherical collapse predictions 
(shaded) capture the trend in the no-chameleon simulations. Pull 
results for n = 1, 2 and consideration of the field evolution suggests 
that spherical collapse predictions should hold for n < 4. Note that 
errors are fully correlated in that the all simulations use the same 
initial conditions and are compared against the same set of ACDM 
simulations. 



this overall enhancement the impact of changing the scal- 
ing parameter n is less significant. This weak dependence 
is in spite of the rapid change in the background fn field 
shown in Fig. [TJ 

We can understand this relative insensitivity by com- 
paring the full simulation results to the no-chameleon 
simulations where the Compton wavelength is fixed to 
its background value through Eq. (|10|). Mass function 
enhancements in the chameleon and no-chameleon sim- 
ulations are nearly the same up until the very highest 
masses. For the large field value today |/ro| = 10~ 4 , 
cluster potential wells are not sufficiently deep to mani- 
fest the chameleon mechanism today. The small effect at 
the very highest masses in fact comes from the chameleon 
mechanism becoming effective at high redshift as we shall 
see. One can in turn understand the relative insensitiv- 
ity to n in the no-chameleon simulations by examining 
the background Compton wavelength evolution in Fig. [T] 
Note that for n < 4, the Compton wavelength varies little 
for rcdshifts z < 1. 

The spherical collapse predictions outlined in the pre- 
vious section are also shown in Fig. [2j The upper bound- 
ary of the shaded region represents enhancements pre- 
dicted by the unmodified spherical collapse parameters 
A v = 391 and S c = 1.673 whereas the lower boundary 
takes the modified parameters A v = 309 and S c = 1.692 
(Table HI]). 

The spherical collapse predictions model the results 
equally well for the n = 1 and n = 2 models. In the high 
mass cluster regime, the unmodified parameters match 
the simulations better. In the low mass end the modi- 
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FIG. 4: Mass function enhancement at z = 0.316 in the large 
field l/ijol = 10~ 4 for n = 1,2. Fractional enhancements at a 
fixed mass remain significant at higher redshift due to the increased 
rarity of such halos in ACDM and the trend remains well captured 
by spherical collapse predictions (shaded region). Field evolution 
in the n = 2 case makes the chameleon suppression in the full 
simulations moderately more important. 



fied parameters agree better. The modified parameters 
also provide conservative estimates of the enhancements 
across the full mass range [IH . 

We further test the large-field no-chameleon simu- 
lations against spherical collapse predictions for even 
steeper n models in Fig. [31 These predictions, based 
mainly on the instantaneous linear growth function, re- 
main accurate despite the extremely strong scaling of the 
force modification with redshift in these models. Further- 
more, Fig. [T] implies that the no-chameleon results should 
be a reasonable approximation to the full simulations for 
n < 4. Thus, in the large field regime, one way to map 
cluster constraints obtained at a given mass M v on one 
/(i?) model to another is to match the linear variance 
cr(M v ). A better approximation can be obtained by set- 
ting the mass function n\ n ^j v equal as we shall see in 

info 

In Fig. 21 we show the mass function enhancements 
at an intermediate redshift z = 0.316 for the large field 
model. Note that the abundance of halos of mass M at 
z = is equal to that of halos of mass M/1.5 at this red- 
shift in a ACDM model, due to the evolution of the mass 
function, and we have adjusted our binning to take this 
into account. Thus for a fixed mass, the enhancement in 
the cluster abundance remains significant. Interestingly, 
the range in spherical collapse predictions continues to 
model these trends once the collapse parameters are ad- 
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FIG. 5: Mass function enhancement in the small field regime. The 
chameleon effect suppresses the enhancement when the background 
field amplitude |/i?ol drops below the depth of the gravitational 
potential for an object of mass M20O' Comparison of the full and 
no-chameleon simulations shows that the limiting mass at which the 
chameleon appears scales roughly as expected: M c h a m °c |/_ro| 3//2 
nearly independently of the scaling index n. Spherical collapse 
predictions roughly capture this suppression in the cluster regime 
-M200 J 3 X 10 14 Mq /h but fail to model the enhancement below 



justed to the matching redshift (see Tab. Hi]). 

The n = 2 results at z = 0.316 show a slight increase 
in the importance of the chameleon suppression when 
compared to z = or n = 1 at the same redshift. This 
is consistent with the suppression of the field amplitude 
shown in Fig. [TJ For n = 2, the effect is only a small 
fractional contribution and spherical collapse predictions 
still work well but suggest that the no-chameleon approx- 
imation may have a smaller range of validity in n at high 
redshift. More generally modified gravity models which 
possess this type of non-linearity that suppresses devi- 
ations in high density regions typically do not predict 
larger enhancements of the cluster abundance at high 
versus low redshift at a fixed degree of rarity or peak 
height v [H]. 



\f R0 \ = 10- 5 , 10" 6 and, for n = 2, \f R0 \ = 3 x lO" 6 . The 
first thing to note is that in the no chameleon simulations 
the impact of changing the field value from \fn,o\ = 10~ 4 
to |/flo| = 10 ~ 5 is less than a factor of 2 in the abun- 
dance at the highest mass bin. We shall see in the next 
section, that this logarithmic sensitivity translates into a 
strong model dependence of observational constraints on 
the field amplitude and Compton wavelength. 

Small field results show a large impact from the 
chameleon suppression as can be seen by comparing the 
full simulations to the no-chameleon simulations. A halo 
is chameleon-screened whenever its gravitational poten- 
tial is larger than the field amplitude in the background 
\fno\- This can be used to derive a threshold mass for 
chameleon screening at z = for a given value of fno 
(see [23j]). We then expect the mass scale M c h a m of the 
chameleon suppression in the mass function to scale sim- 
ilarly as the threshold for chameleon screening. In par- 
ticular, M c ham should depend mainly on |/ro| and only 
weakly on the scaling index n. Specifically, neglecting 
the small mass-dependence of the halo concentration, we 
would expect the onset of the chameleon suppression to 
scale as M cham cx \fno\ 3/2 - 

We see from Fig. [5] that the results are consistent with 
this scaling: roughly, the chameleon for |/i?o| = 10~ 5 is 
important for A/200 > 6 x lO 14 M //i while for \/rq\ = 
10~ 6 , the suppression appears at M200 > 2 x lO 13 M //i. 
The I /ho I = 3 x 10~ 6 , n = 2 case falls consistently right 
in between the two despite being a different n model. 

Spherical collapse predictions roughly model the re- 
duced enhancement in the cluster regime of M200 > 3 x 
1O 14 M jh. They correctly predict an absence of a signifi- 
cant enhancement for \jm\ ^5 3x 10 -6 . However, unmod- 
ified collapse parameter predictions can fractionally over- 
estimate the enhancement unlike in the large field regime, 
while modified collapse parameter predictions predict a 
reduction in the cluster abundance (Ani n m < 0) not seen 
in the simulations. Moreover both cases do not predict 
the correct behavior at lower masses where the full sim- 
ulations possess a higher abundance of halos than both 
the no-chameleon simulations and the collapse predic- 
tions. Hence in the small field regime they should not be 
used for constraints from galaxy groups or smaller mass 
objects or if precision predictions are required at cluster 
masses. We defer such modeling to a future work. 



Current Constraints 



B. Small Field Regime 

As cluster abundance and other cosmological tests im- 
prove, the large-field models will be excluded (if no order 
unity excesses over ACDM expectations are detected). 
In the small field regime of |/ro| < 10 ~ 5 , the chameleon 
mechanism is effective even today. 

In Fig. [5J we show small field results for n — 1 and 



Given that spherical collapse predictions work equally 
well for all of our broken power law models with n < 4 
in the cluster regime and capture the scalings seen in the 
full simulations, we can extend the constraints on the 
f(R) model with n = 1 Q that were obtained using the 
observed abundance of local X-ray clusters selected in the 
ROSAT All-Sky Survey and followed up with Chandra 
observations 24 1. 

The constraints were obtained by using the spherical 



7 



0.1 H I I I | I I I I | I I I I | I I I I | I I I I | I I H 




n 



FIG. 6: Constraints on |/i?o| (upper panel) and the Compton 
wavelength Ac (lower panel) as function of the index n. We 
have converted the 95% confidence level upper limits on |/ao 
reported in [l[ for n = 1 to other values of n using the spher- 
ical collapse model as described in the text. The medium 
shaded band corresponds to the default limit reported in [lj|, 
while dark and light shaded areas use more or less conserva- 
tive assumptions (see text). 



collapse model (see section !!! C|) to predict the f(R) mass 
function enhancement at a pivot mass of Mx,eS ~ 3.7 x 
10 14 MqIK 1 for an overdensity of 500 with respect to 
critical density. Fig. [3] shows that the spherical collapse 
model is equally valid for other values of n as long as the 
chameleon effect is negligible, and it is straightforward to 
translate the constraints to other values of n by matching 
the abundance at Mx, c fi- 

The results are shown as function of n in Fig. [5] for 
a range of conservative to aggressive interpretations of 
the data and modeling (see [l| for further discussion). 
In the top panel we show the 95% statistical limits on 
the field amplitude today fuo and in the bottom panel 
the Compton wavelength in the background today Xco- 
The medium shaded region shows the result for the de- 
fault constraint, \fno\ < 1-3 x 10~ 4 at n = 1, using 
the modified spherical collapse parameters (lower edge 
of shaded band in Fig. The dark region shows the 
most conservative constraints (|/ito| < 3 x 10 ~ 4 ), us- 
ing the modified collapse parameters and in addition 
assuming X-ray masses are underestimated by 9%. Fi- 
nally, the light region shows more aggressive constraints 
(l/flol < 4 x 10~ 5 ), using the unmodified collapse param- 
eters (upper edge of shaded band in Fig. [3]). Note that 
even this case is still somewhat conservative, since for 
clusters at fixed mass, dynamical mass estimates such as 



X-ray masses will be enhanced by ~ 20% in the large-field 
limit of f(R) gravity HH, due to the increased depth of 
the potential well. This increases the abundance at fixed 
Mx in f{R) considerably. 

While the change in the fractional enhancement of the 
mass function from ACDM with n is relatively small, the 
impact on the model parameters can be large. Specif- 
ically between the n = 1 and n — 4 models the field 
amplitude limits change by over an order of magnitude 
and Compton wavelength constraints by a factor of sev- 
eral. 

Nonetheless the cluster abundance measurements can 
already rule out a substantial portion of the cosmolog- 
ically interesting regime for all cases, limiting the al- 
lowed range of enhanced forces to 10 — 100 Mpc. Fu- 
ture large cluster samples have the potential to push the 
limits down by an order of magnitude before chameleon 
effects cause a suppression of the enhancement. 



IV. DISCUSSION 

We have conducted TV-body simulations to test the en- 
hancement of the cluster abundance in a variety of f(R) 
models. These models differ in the redshift evolution 
of both the linear force enhancement and the non-linear 
chameleon mechanism which suppresses such enhance- 
ments in the deep gravitational potential wells of clusters 
of galaxies. These results test the robustness of model in- 
dependent techniques such as spherical collapse for pre- 
dicting the enhancement and constraining modified grav- 
ity with cosmological data. 

We find that for cluster mass halos, the spherical col- 
lapse predictions work equally well for different models 
at least as long as the redshift evolution of the field is not 
so steep as to invalidate the division between large field 
and small field regime imposed at \fm\ ~ 10~ 5 for the 
background field amplitude at z = 0. In the large field 
regime the background field amplitude is larger than the 
depth of the gravitational potential wells of clusters and 
hence the chameleon effect is inoperative. For a scaling 
index n < 4, a large field model retains this property 
for z < 1 when clusters form. In this regime, the frac- 
tional enhancement of the cluster abundance relative to 
ACDM is a relatively weak function of n that is deter- 
mined by the evolution of the Compton wavelength or 
range of the force in the background. In the opposite 
small field regime, the enhancements become suppressed 
above a limiting mass that depends mainly on the field 
amplitude M cha m cx |/i?o| 3/2 - 

We use these results to extend the implications of the 
local cluster abundance to the whole class of broken 
power law models. Most models in the literature can 
be mapped onto the parameter space of this class. Con- 
straints on the field amplitude and Compton wavelength 
today are strongly model dependent due to the logarith- 
mic dependence of the cluster abundance on their values 
in any given model. 
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Results based on different model assumptions can be 
mapped onto each other by matching instead the lin- 
ear theory rms fluctuation at the radius implied by the 
observed mass scale or even more directly by matching 
spherical collapse mass function predictions as we have 
shown for the local X-ray cluster abundance. 
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